library(logr)
lf<-log_open("03_simulation_appendix.log")
#---------------------------------------------------------------------------------------------------
#APPENDIX
#---------------------------------------------------------------------------------------------------
#TABLE A-3.1:
df_coverage_se= results_all %>%
    group_by(blocking, n_sample) %>%
    filter(blocking != "Scenario 3: Block on Outcome-Related Variables") %>%
    summarize(
        iv_se = mean(iv_se),
        iv_cov = mean(CACE > (iv- 1.96 *iv_se) & CACE < (iv + 1.96 * iv_se)),
        iv_block_se = mean(iv_block_se),
        iv_block_cov = mean(CACE > (iv_block - 1.96 *iv_block_se) & CACE < (iv_block + 1.96 * iv_block_se)),
        block_dim_se = mean(block_dim_se),
        block_dim_cov = mean(CACE > (block_dim - 1.96 * block_dim_se) & CACE < (block_dim + 1.96 * block_dim_se)),
        psw_se = mean(psw_se),
        psw_cov = mean(CACE > (psw - 1.96 *psw_se) & CACE < (psw + 1.96 * psw_se)),
        block_psw_se = mean(psw_block_se),
        block_psw_cov = mean(CACE > (psw_block - 1.96 * psw_block_se) &
            CACE < (psw_block + 1.96 * psw_block_se)),
        placebo_se = mean(placebo_srs_se),
        placebo_cov = mean(CACE > (placebo_srs - 1.96 * placebo_srs_se) &
            CACE < (placebo_srs + 1.96 * placebo_srs_se)),
        block_placebo_se = mean(placebo_block_se),
        block_placebo_cov = mean(CACE > (placebo_block - 1.96 * placebo_block_se) &
            CACE < (placebo_block + 1.96 * placebo_block_se)),
    )
log_print(xtable::xtable(df_coverage_se[,-1]), include.rownames=FALSE)


#TABLE A-3.2:
log_print(xtable::xtable(data.frame(df_mse[,-1], df_bias[,-1])), include.rownames=FALSE)

#------------------------------------------------------------------------------------------
#ER VIOLATION:
load('Simulation Results/sims_er_violation.Rdata')
df_mse_er <- results_er %>%
            group_by(block_var, beta) %>%
                summarize(
                    #at = mean((AT - CACE)^2),
                    iv=mean((iv-CACE)^2),
                    iv_block=mean((iv_block - CACE)^2),
                    placebo = mean((placebo_srs - CACE)^2),
                    placebo_block = mean((placebo_block - CACE)^2),
                    psw = mean((psw-CACE)^2),
                    psw_block = mean((psw_block - CACE)^2),
                    bdim = mean((block_dim - CACE)^2)
                )


df_var_er <- results_er %>%
            group_by(block_var, beta) %>%
                summarize(
                    #at = mean((AT - CACE)^2),
                    iv=var(iv),
                    iv_block=var(iv_block),
                    placebo = var(placebo_srs),
                    placebo_block = var(placebo_block),
                    psw = var(psw),
                    psw_block = var(psw_block),
                    bdim = var(block_dim)
                )
df_plot_er = rbind(
    data.frame(df_mse_er %>% gather(key = Estimator, value = Value, -c(block_var, beta)), Measure="Bias^2"),
    data.frame(df_var_er %>% gather(key = Estimator, value = Value, -c(block_var, beta)), Measure="Variance")
    )

df_plot_er = recode(df_plot_er)

#FIGURE A-3.2 (Estimator performance under Exclusion Restriction Violation)
df_plot_er %>%
filter(!Estimator %in% c("Placebo (CR)", "Placebo (Block)")) %>%
ggplot(aes(x = as.factor(beta), y =Value, fill=Estimator, alpha=Measure, group=Estimator))+
geom_col(position='dodge')+
facet_wrap(~block_var)+
scale_alpha_manual(values = c(0.5, 1))+
ggtitle("MSE under Violations of Exclusion Restriction")+xlab("Magnitude of ER Violation")+
scale_fill_manual(values=pl1)+theme(axis.text=element_text(size=12)) + ylab("Mean Squared Error")
ggsave("../Figures/figA3_2_mse_breakdown_er.pdf", width=11, height=4)

log_close()